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We present a discrete analysis of non-reflecting boundary conditions for the discon- 
tinuous Galerkin method. The boundary conditions considered in this paper include the 
recently proposed Perfectly Matched Layer absorbing boundary condition for the linearized 
Euler equation and two non-reflecting boundary conditions based on the characteristic de- 
composition of the flux on the boundary. The analyses for the three boundary conditions 
are carried out in a unified way. In each case, eigensolutions of the discrete system are ob- 
tained and applied to compute the numerical reflection coefficients of a specified out-going 
wave. The dependencies of the reflections at the boundary on the out-going wave angle 
and frequency as well as the mesh sizes are studied. Comparisons with direct numerical 
simulation results are also presented. 


Introduction 

The discontinuous Galerkin method (DGM) provides an 
attractive platform for computational aeroacoustics in 
dealing with complex geometries while using high or- 
der approximations [3j. Our recent works have shown 
that discontinuous Galerkin scheme is super-accurate 
for numerical simulation of wave propagations [6,7]. A 
framework of discrete analysis has been proposed in [6,7] 
to study systematically and analytically the numerical 
errors occurred due to a sudden discontinuity in mesh 
topology or order of the basis functions in one and two 
spatial dimensions. In this paper, this framework is ex- 
tended to study numerical errors due to non-reflecting 
boundary conditions in DGM schemes. The effects of 
various numerical non-reflecting boundary conditions 
will be treated in a unified way as either a variation 
in the flux formula or a variation in the underlying gov- 
erning equations across a mesh interface. 

Non-reflecting boundary condition is an essential part of 
any numerical code in computational aeroacoustics. In 
this paper, we will analyze three types of non-reflecting 
boundary conditions for DGM, namely, the characteris- 
tic non-reflecting boundary condition, finite wave model 
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boundary condition and the Perfectly Matched Layer 
(PML) absorbing boundary condition. An earlier ver- 
sion of PML for the linearized Euler equation was given 
in split physical variables. Its implementation in the 
discontinuous Galerkin schemes has been demonstrated 
in [1], An unsplit version has been proposed recently in 
[5] that eliminates numerical instabilities that can oc- 
cur in the earlier version. The implementation and the 
discrete analysis of the unsplit version will be studied in 
the present paper. We will recognize the advantage of 
DGM in that, unlike the finite difference schemes, the 
absorption coefficient in the PML domain is no longer 
required to be continuously varying but can be increased 
discontinuously. Although the PML is reflectionless 
in its formulation in the non-discrete partial differen- 
tial equations, the discretization of the equations can 
nonetheless cause reflections. We will first present an 
eigensolution analysis of the discretized PML equations 
in DGM. Then, the results of the eigensolution analysis 
will be applied to study wave reflections caused by the 
discretization process at the interface of the Euler and 
PML domains. 

We will also present a discrete analysis of the character- 
istic non-reflecting boundary condition [8,9] as well as 
the finite wave model boundary condition proposed in 
[2]. The characteristic non-reflecting boundary condi- 
tion is widely used in DGM schemes because it is easy 
to implement due to the intrinsic upwinding features 
of the schemes. The finite wave model non-reflecting 
boundary condition proposed in [2] improves the char- 
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acteristics boundary condition in that the reflection of 
out-going waves can be annihilated under certain con- 
ditions. The results of our discrete analysis will be 
compared with the theoretical reflection coefficients ob- 
tained from the non-discrete partial differential equa- 
tions as well as those obtained from direct numerical 
simulations. 

In the next section, we will first describe the discretiza- 
tion of the PML equation in the DGM schemes and 
present an eigensoultion analyses of the semi-discrete 
system. Since the Euler equation is a special case of 
the PML equation, the discrete analysis can be readily 
extended to the analysis of the characteristic and finite 
wave model non-reflecting boundary conditions in sub- 
sequent sections. 

Discrete analysis of PML absorbing 
boundary condition 

It has been shown in [5] that the PML equation is per- 
fectly matched to the Euler equation (1). That is, in 
the non-discrete form, there is no reflection created at 
an interface between the Euler and PML domains. The 
effects of numerical discretization by the DGM scheme 
on this reflectionless property will be studied in two 
steps. First, we will compute the eigensolutions of the 
semi-discrete PML equation and study the effects of the 
absorbing coefficients oh the wave propagation proper- 
ties. Second, we will utilize the eigensolutions so formed 
to compute the wave reflection and transmission coef- 
ficients at an interface of discretized Euler and PML 
domains. 

We start with the the linearized Euler equation in con- 
servation form, 


dxx 

dt 


+ V • F(u) = 0, 


where F(u) is the flux vector given by 


(1) 


xi and x -2 directions respectively; p is the pressure; M 
is the Mach number of the mean flow. For conciseness, 
the density equation is not included here. 

At non-reflecting boundaries, absorbing layers are uti- 
lized so that the waves exiting the physical domain are 
damped. The PML equations to be used in the absorb- 
ing domain can be written as follows [5], 


— + V • F pm ;(u, q) + (a x + cr y )u + a x a y q 


+cr x PA 1 (u + a v q) =0, (4) 


and 


<9q 

dt 


= u, 


(5) 


where q is a vector of auxiliary variables and /3 in (4) 
is a function of the mean flow Mach number as 



M 

1 - M 2 ' 


The flux vector F pm / in (4) is of the form 

Fpmi(u,q) = (Ai (u + er y q), A 2 (u + ff*q)). (6) 

Here the absorption coefficients a x and <r y can be con- 
stants or functions of x and y respectively [5]. 

We consider a discontinuous Galerkin scheme for (4)-(5) 
in which the spatial domain is partitioned into elements 
denoted by fl nm , where n and m are the element indices 
in the directions of xi and x 2 (spatial coordinates) re- 
spectively, as in Figure 1. For simplicity, we will only 
consider uniform rectangular elements in this paper. 
Triangular elements can be studied in a similar ap- 
proach as we have done in [7] for the Euler equations. 
For each element in the PML domain, the numerical 
solution vectors, denoted by u™ n and q£ m , are approx- 
imated by an expansions in polynomials 


in which 


F(u) = (Aju, A 2 u), 


(2) 


'(x,f)=E c rwrw, 


(7) 


f=0 
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/ u \ 
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0 0 0 
a 2 = j 0 0 1 
0 1 0 


(3) 


where u and v are velocity components (non- 
dimensionalized by the speed of sound) in the spatial 


q r 


(x,«) = £drwrw, 


(8) 


<=0 


where {0J? m (x),f? = 0, 1, .., L} is the set of the basis 
polynomials and c" m (f) and d" m (t) are the expansion 
coefficients for u£ m and q£ m respectively. A weak for- 
mulation of (4)-(5) is obtained by substituting (7)-(8) 
into (4)-(5) and requiring that the equations be orthog- 
onal to the basis functions. After applying integration 
by parts, the semi-discrete equations for (4) and (5) are 
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f ft ii nm f 

/ lh^(x) rf x + / (nF pml )<tf m (x)ds 

*^n n m ^r nm 


-/ 


and 


F pm j • V4>r(x)dx + f s r^? m (x)dx = 0, (9) 

J ttnm 



dq r 
dt 


<f>? m (x)dx 



Uh m </>? m (x)dx, 


( 10 ) 


for f? = 0, 1,...,L, where r nm denotes the boundary of 
fl nm and n = ( ni,ri 2 ) is its outward normal. Here s£ m 
in (9) represents the non-differential terms in (4): 

S r = (a x +ay)ur+v x cj y qr+v x 0A 1 (ur+°y< m )- 


Equations (9) and (10) will result in a system of ordi- 
nary differential equations for the expansion coefficients 
c" m (t) and d" m (f). 

For the PML flux vector given in (6), the normal flux 
at the edge of an element is 


n • F pm i = (ri! A! + n 2 A 2 )u + (cr y niAi + a x n 2 A 2 )q. 

( 11 ) 

For convenience of discussion, we denote the normal flux 
in the original Euler equation by ' ' 

n •. F = (nj Ai 4- n 2 A 2 )u = A n u, (12) 

which is the same as the first term in n • F pm ; above. 
Further, it is easy to verify that the second term in (11) 
can be simplified as 


evaluated using the solutions inside and outside of that 
element respectively. Then both flux formulas can be 
conveniently written as follows 

n • F pmi = A + (u+ + aq+ ) + A_ (u^ + crq^ ) (15) 

where 

A+ = i(A n + |A„|) and A_ = ^(A n - |A n |) (16) 

for the characteristics splitting flux formula and 

A + = ~(A n -|- |Umaj:|I) Und A„ = — ( A n |Omax|I) 

(17) 

for the Lax-Friedrichs flux formula, where a max is the 
largest (in magnitude) eigenvalue of A n . 

In our analysis, we assume that the same flux formula 
is used through out the computational domain. 
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( <7;, n i A i + <7x 712 A 2 ) q 


J a y A n q vertical edges (n 2 = 0) 

) a x A„q horizontal edges (r»i = 0) 

(13) 


This leads to the following more compact flux expres- 
sion on the edges 


Figure 1. A schematic of partitioning of the computational 
domain in quadrilateral elements. Here, n and m are the 
indices of the element and the number inside the parentheses 
indicates the ordering of the element edges. 


n-F pmi = A n (u + <rq) (14) 

where a is either a x or a y according to (13). 

Since we do not require that the numerical solution be 
continuous across the interface of any two elements, the 
normal flux n-F pm ; appearing in the edge integral in (9) 
is not uniquely determined. To complete the formula- 
tion, a numerical flux formula needs to be specified. As 
in references [6] and [7], we will consider two commonly 
used flux formulas, namely, the characteristics splitting 
flux formula and the Lax-Friedrich flux formula. For 
convenience of discussion, let u£,q^ and u^,q^ de- 
note the values of u/ ( , q^ on the edge of an element 


We order and denote the four edges of fi nm as rl'i, i = 
1,2, 3, 4, as shown in Figure 1. Using (15) for the normal 
flux, the weak formulation (9) can now be written in 
detail as 

r rh\ nm 

/ fc(x)-^-dx 

+ f <Mx)[ Afur + A(! ) u” m-1 ]ds 

+ [ M^)[A^u n y m +A^n^ +lm ]ds 

Jrlfl 
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+ / <Mx)[ A®ur + AL 3) ur +1 ]ds 

n m 

+ f fc(x)[ A^u"™ + AL 4 ) <- lm ]d s 

-f tp-Ax u; m dx- [ ^p-A 2 u2 m dK 
Jn nm dx Ju nm d y 

+ f sr0? m (x)dx = O (18) 

where we have used notation 

, nm _ nm . nm / in \ 

— U/, + &x,yQh • (19) 

The superscripts (1) — (4) in equation (18) indicate the 
specific edge of the element for the integrals as shown 
in Figure 1. 

Eigensolutions of the discrete system 


u£ m (f, V, t) = e~ iut (P T ® I)C nm , (22) 

<lh m (Z,V,t) =e- iwt (P T ®I)t) nm . (23) 

By substituting (20) and (21), or their tensor forms (22) 
and (23), into the semi-discrete PML equations (18) and 
(10), we get a linear system for the expansion coeffi- 
cients C nm and D nm as follows: 

-iw { Q ® I)C nm + (BP) ® ASICS'" + (B'O) ® A^ 1 ) )C I nm_1 
+(B< 2 > ® + (B'< 2 > <g> A ( _ 2) )C£ +lm 

+ (B (3) ® A^ ) )C” m + (B'< 3 > ® AL 3) )C x nm+1 
+ (B< 4 > ® A< 4 »)c7 m + (B'< 4 > ® AL 4) )c !/ n - lm 
-(Q* ® Ai)C, nm - (Q„ ® A 2 )C x nm 

+ (Q ® I)[(<Tx + cq,)C nm + a x a y t> nm ] + a x P( Q ® Ai)C? m = 0 

(24) 

and 


To look for wave-like solutions to (18), we assume peri- 
odicity in time with a frequency to and let 


L 



(20) 

(=0 


L 


e=o 

(21) 


where £ and rj are local coordinates. The coefficients 
c? m and qj* m are now independent of t. 

For convenience of discussion, define vectors 


-iui(Q ® I)D nTn = (Q ® I)C nm (25) 

where 

C ~ nm f^nm , _ nnm 

x,y — ^ + Vx,ylJ 

and Q is the “mass matrix”, Q x and Q y are the “stiff- 
ness” matrices, and BW and B'b) are the “edge” ma- 
trices resulted from the edge integrals in (18) involving 
solutions inside and outside of the element respectively. 
The detailed definitions of these matrices are given in 
the Appendix (Al). 

By (25), we easily get 


and 



/ c" m \ 


/ d” m \ 

£yrim 

c" m 

,D nm = 

dr 


l c r ) 


\ dr 


i 


\ 


P(^) = 


fcitv) 

V 4 >l(Z,v) j 


D nm = -C nm . (26) 

to 

As a result, equation (24) can be written in C nm only 
as 

-iu{ 1 + — )(1 + — )(Q ® I)C nm + H (0) C nm 

U LJ 

+(i + — )H (l) c nm-1 + (l + ^-)H (2) c n+lm 

U) LU 

+(1 + ^A)H (3) C nm+1 + (1 + ^) H (4) C"- lm = 0. (27) 

U) U) 


Then, using the Kronecker product ® (see reference [6] 
for definition and its properties), the solution in (20) 
and (21) can be expressed as follows, 


The matrices can be obtained by directly com- 

paring (27) with (24) and will not be given explicitly 
here. 
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For a uniform distribution of elements, the coefficient 
matrices in (27) are independent of the element indices 
n and m. We seek solutions of the form 

C nm = A"A™C (28) 

where Aj and A 2 are related to wavenumbers k\ and fc 2 
in the x\ and x 2 directions, respectively, as follows 


Ai = e iklhl , \ 2 = e ik ^ (29) 

where h\ and ft 2 are the spatial dimensions of elements 
in £1 and £2 directions respectively. An eigenvalue 
problem for C is formed by substituting (28) into (27) 
which yields immediately 


-*w( 1 + — )(1 + — )(Q ® I)C + H(°)C 

UJ U) 

+(1 + — + Ai(1 + ^)H (2) C 

U) A 2 U) 

+(1 + — )A 2 H (3) C + -1(1 + = 0. (30) 

oj Ai u) 

To study the spatially propagating waves, we solve for 
Ai and C as the eigenvalue and eigenvector from (30) 
for given a/, the frequency, and fc 2 , the wavenumber in 
£2 direction. The corresponding discrete wave number 
ft* of the eigensolution can be found according to (29) : 


Jfeffti = —i ln(Ai). 


It is well known that, for any given values of w and 
k 2 , the Euler equation (1) supports two acoustic waves 
given by the dispersion relation 

(w - Mki) 2 - k\ - ftf = 0 (31) 

and one vorticity wave given by 


u - Mki = 0. (32) 

The acoustic waves can be conveniently expressed in the 
given frequency ui and a wave angle <p as 


w cos 4> oj sin (j> 

1 + M cos <p ’ 2 1 + M cos <f> 


(33) 


To illustrate typical eigensolutions of (30), we solve the 
eigenvalue problem with square elements h\ = h 2 = ft 
and the order of basis functions p — 2. In Figure 2, we 
show the discrete wavenumber k[ corresponding to the 
right-traveling acoustic mode obtained as the eigenvalue 
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of (30). Plotted in Figure 2a are the ratio of the real 
part of ft* and the exact fti in (33) as a function of non- 
dimensionalized frequency uih. In Figure 2b, we show 
the imaginary part of ft] 1 . The value of ft 2 , hence the 
value of A 2 in (30), is given according to (33) and (29) 
with <j> = 7r/6. The Mach number of the flow is M = 0.5. 
For simplicity, a y = 0 in all the calculations. The solid 
line in both figures represents the result for a x — 0, 
i.e., that of the discretized Euler equation, and dashed 
lines represent that of the PML equation with non-zero 
a x as indicated. Figure 2 illustrates the effects of the 
absorption coefficient on the discrete wavenumber in the 
PML domain. Figure 2a indicates that an increase in 
the absorption coefficient has very little effect on the 
real part of ft* for waves within the resolved range in 
non-dimensionalized frequency u>h < 1, see reference 
[7]. The imaginary part of ftj, however, has an increase 
that is proportional to the increase in the absorption 
coefficient, indicating exponentional absorption of the 
wave. Furthermore, the damping rate remains the same 
for all frequencies within the resolved range cuft < 1. 
This means that the absorption rate for the acoustic 
wave will be independent of the wave frequency. 

In Figure 3, we show the dependence of PML damping 
rate as a function of the wave angle for a given wave 
frequency tuft = 0.5. In general, the damping rate re- 
duces as the wave angle 6 increases. This is consistent 
with the theoretical prediction of the non-discrete PML 
equation [5]. 

Aside from the physical modes, there are also non- 
physical eigensolutions, or spurious modes, for (30). 
The number of spurious modes depends on the order 
of the basis functions and the flux formula being used. 
It can be shown that there are p spurious modes for 
each physical wave of the non-discrete partial differen- 
tial equation when the exact characteristics splitting is 
used in the flux formula. The number of the spurious 
mode will be increased to be 2p + 1 when Lax-Friedrichs 
flux formula is used. The non-physical waves are irreg- 
ular and highly damped [6,7]. 

Reflection and transmission at an interface of Euler 
and PML domains 

Although the PML is reflectionless in the non-discrete 
limit, discretization of the partial differential equations 
may cause numerical reflections at the interface of the 
Euler and PML domains. We will now analyze the nu- 
merical reflection error due to the discretization process. 

Without loss of generality, suppose that the elements in 
the Euler domain have indices n < 0 and those in the 
PML domain have indices n > 1, as shown in Figure 
4. We will introduce an out-going wave (right-traveling 
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Figure 2. Real and imaginary parts of the numerical 
wavenumber k' of the physical mode of (30) as functions 
of wave frequency wh. h\ = h% ~ h, M = 0.5, p = 2. 
u y = 0. 



wave angle <f> 

Figure 3. Normalized numerical wave speed and spatial 
damping rate of the physical mode as functions of wave 
frequency. 9 = 7r/6, p = 2, h\ = = 1. o v = 0. 

uih = 0.5. 
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Figure f. A schematic of wave reflection and transmis- 
sion at an interface of Euler and PML domains. The 
<j>i is the wave angle of the out-going wave. 


in Figure 4) in the Euler domain and compute the re- 
flected and transmitted waves in the Euler and PML 
domains respectively. For an interface shown in Figure 
4, the reflected waves will include all the left-traveling 
eigensolutions of the discrete Euler equation and the 
transmitted waves will include all the right-traveling 
eigensolutions of the discrete PML equation. Both the 
physical and non-physical modes have to be considered. 
Specifically, we assume the solution in each half domain 
as 

for n < 0: 

ClZer = vr« + PwlVWI + E 0B* , (34) 

3 

for n > 1: 

C ~ nm apml^Tnm , \ ' opml^rnm /OK - * 

pml ~ Pw « V W R + Pe r V E r > (35) 

3 

where subscripts W R and W L denote the physical right 
and left traveling waves respectively, and E R and Ej 1 
denote the non-physical right and left traveling spurious 
wave modes respectively. Here the V vectors are the 
corresponding eigenvectors of the semi-discrete equa- 
tions. In (34) and (35), the out-going wave, V^J,, has 
an amplitude of unity, and (3 w l and /3 e l are the reflec- 
tion coefficients and /?vy/?and (3 e r are the transmission 

j 

coefficients. 

The coupling of the two solutions occurs at elements 
immediately on either side of the interface, namely, the 
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elements with indices n = 0 and n = 1. The semi- 
discrete equation for the Euler domain is that given in 
(27) with a x = a y = 0. For the PML domain shown in 
Figure 4, we have a y — 0 but a x 0. Thus, we have 

for n — 0: 


— icu(Q ® I)C' 


0,m 

Euler 


4- H 


(o) 


, it 

^ Euler + “ 


( 1 ) 


0 , m — 1 
Euler 


+H 


(2) 


C. + H (3) l 


iO,m+l 

'Euler 


+ H 


(4) 


f-i — l,m 

^ Eui 


= o, 


(36) 


for n = 1: 


— Mi + -f)(Q®i)cir, + H w c. 


r (0)^l,m 
''pml 


+( i +7 )h(1) ^t'+ h(2) c;: 

+(1 + — )H (3) C; m r l + H (4) C^ er = 0. (37) 

UJ 

Here, the absorption coefficient a x is assumed to be 
constant for all the elements in the PML domain. 


Equations (36) and (37) can be further simplified by the 
fact that C|” l ifr and are formed by the eigenvec- 
tors. As. a result, we get the following more compact 
matching conditions [7], 

H (2) C^ er = H( 2 )C^, ' (38) 

■ H (4) C° £ l=H( 4 )c; m m ,. (39) 


By substituting (34) and (35) into (38) and (39), we can 
find the numerical reflection coefficients at the interface. 

As an example, we consider an out-going acoustic wave 
with a wave angle <j>i = 7r/6. The order of the basis func- 
tions used is p = 2. For this order, when the exact char- 
acteristic splitting formula (16) is used, the reflected 
waves in the Euler domain include the left-traveling 
acoustic mode, denoted by W L , and two non-physical 
spurious modes, denoted by Ef and Ef. Figure 5 shows 
the reflection coefficients of the physical {W L , in cir- 
cles) and non-physical (Ef , E.f, in diamonds) modes 
as functions of the non-dimensional wave frequency uih. 
As the value of uih decreases, i.e., when the frequency of 
the wave is reduced (hence its wavelength is increased) 
or the size of the element is decreased, the reflection at 
the interface caused by numerical discretization also de- 
creases. More significantly, the reflection coefficient for 
the physical mode is decreasing at a higher order than 
that of the non-physical modes. Since the non-physical 
modes are highly damped by the numerical dissipation, 
Figure 5 indicates that the reflections occurred at the 
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Figure 5. Reflection coefficients of the physical (W L ) 
and non-physical (Ef , Ef) modes with basis function 
order p = 2. M = 0.5, o x = 1. s indicates slope, of 
solid line. The out-going acoustic wave angle <fi = 7r/6. 
hi = /12 — h. 



Figure 6. Reflection coefficient of the physical mode as 
a function of wave angle. The width of the PML domain 
is D . p = 2. M = 0.5, G x h = 1. u>h = 0.5. 


interface will likely have only a local effect and will not 
be propagated globally to the computational domain. 
We note that when the Lax-Friedrichs flux formula was 
used, the number of reflected non-physical modes was 
increased. However, trends similar to those of Figure 5 
were obtained. 

The transmitted waves are damped exponentially in- 
side the PML domain, at a rate that is determined by 
the imaginary part of the numerical wavenumber k[ . 
In practical computations, however, the PML domain 
is of finite width and wave reflection can occur when 
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the PML domain is truncated. The waves reflected 
at the end of the PML domain will be damped again 
when traveling back through the PML domain before 
re-entering the Euler domain. The magnitude of the 
wave that re-enters the Euler domain relative to that of 
the initial out-going wave will be referred to the total 
reflection coefficient of the PML domain. Obviously, 
the total reflection coefficient depends on the width of 
the PML domain as well as the value of the absorption 
coefficient used. In Figure 6, we show the total reflec- 
tion coefficient of the physical mode as a function of 
the out-going wave angle when the width of the PML 
domain is 2, 4 and 6 elements. The value of the absorp- 
tion coefficient a x h is unity. Here, we have assumed 
complete reflection at the end of the PML domain. The 
reflection can be further reduced, in fact, by using other 
type of non-reflecting boundary condition to terminate 
the PML domain itself. 


Characteristic Non-reflecting Boundary 
Condition 



Figure 7. Schematic of elements at the non-reflecting 
boundary. The elements at the boundary have index n = 
0. (l)-(4) indicates the order of the edges of element 
flnm • 


We will now study the characteristic non-reflecting 
boundary condition applied to the discontinuous 
Galerkin schemes. The implementation of the boundary 
condition is relatively simple by using characteristics 


splitting flux formula (16) at the boundary elements. 
The non-reflecting boundary condition is imposed by 
eliminating the exterior flux contribution in the edge 
integrals along the non-reflecting boundary [8,9]. 

Consider a non-reflecting boundary shown in Figure 7, 
where the boundary elements have index n = 0. The 
semi-discrete equation for the Euler domain is the same 
as those given in (18) or (24) for the PML domain 
with o x = o y = 0. For elements at the non-reflecting 
boundary with index n = 0, the edge integral along the 
non-reflecting boundary (i.e., side (2)) in equation (18), 

f m <M x )[AV 2 ) u°’ m + A ( i ) u^ ,m ]ds, 

is replaced by 

<Mx)A' 2) u°’ m ds 

where A+^ is the flux matrix from the characteristics 
splitting formula given in (16), i.e., 

Af = i(A< n 2 )-HA(% 

We note that when the Lax -Friedrichs flux formula is 
used in the interior elements, Ay ; differs from Ay that 
is used in the eigensolution analysis of (24). They are 
of course the same when the characteristics splitting 
flux formula is used in the interior domain. The total 
reflection error can be further reduced by increasing the 
width of the PML domain. 

In terms of the solution expansion coefficient vector 
C nm , the characteristic non-reflecting boundary condi- 
tion leads to the following equation for the elements at 
the boundary, 

(B ( 2 ) ®A ( f 2 , )C 0 ' m + (B' ( 2 ) igiA) 2 ) )C 1 ’ m = (B ( 2 ) ®Ay 2 ) )C°' m . 

(40) 

By considering a numerical solution of the form similar 
to (34) for C n,m in (40) that consists of an out-going 
physical wave mode and reflected acoustic and other 
non-physical modes, we can compute the reflection co- 
efficients using (40). In Figure 8, we show the magni- 
tudes of reflected physical (acoustic) and non-physical 
wave modes in response to an out-going (right-traveling) 
acoustic wave that has the angle <pi = 7r/6. The order of 
the basis functions used is p = 2 and the characteristics 
splitting flux formula is used through out the compu- 
tational domain. As the wave frequency oj or the mesh 
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Figure 8. Reflection coefficient of the characteristics 
boundary condition. W L is the physical mode (circles) 
and E^and E% are non-physical modes (diamonds). 
M = 0.5. Dashed line is the prediction according to 
the partial differential equation, s indicates the slope of 
the solid line. 


size h decreases, the magnitudes of the reflected non- 
physical modes (in diamond symbols) are reduced and 
approach to zero in the limit. However, the magnitude 
of the reflected acoustic mode (in circles) is approaching 
a non-zero limit predicted by the non-discrete partial 
differential equation model given in the Appendix. This 
is in contrast to the case of PML boundary condition 
shown in Figure 5 where the reflected physical mode 
goes to zero in the limit. It indicates that the reflection 
error can not be eliminated by mesh refinements for the 
characteristic non-reflecting boundary condition. Since 
a well-resolved physical mode has very small numerical 
damping [7], the reflection error can have a global effect 
on the numerical solution. 

In Figure 9, we show the magnitude of the reflected 
physical mode (left-traveling) relative to an out-going 
(right-traveling) acoustic wave. Plotted are the reflec- 
tion coefficients as functions of fa, the out-going wave 
angle, for the mean flow Mach number M = 0 and 0.5. 
At the resolved frequency wh = 0.5, the numerical re- 
flection coefficient coincide with that predicted by the 
non-discrete model shown in dashed lines. Clearly, the 
characteristic boundary condition is most effective when 
the wave angle is close to zero. 

Finite Wave Model 

As we have seen in the previous section, the characteris- 
tic non-reflecting boundary condition works best when 
the out-going wave angle </>, is zero or close to zero, 



out-going acoustic wave angle <t>, (degrees) 

Figure 9. Reflection coefficient of the characteristics 
boundary condition. Circle: M = 0; square: M = 0.5. 
l oh = 0.5. Open symbols are the results of discrete anal- 
ysis and closed symbols are results from direct numerical 
simulation. Dashed line is the prediction according to 
the solution of partial differential equations. 


i.e., when the direction of wave propagation is nearly 
normal to the boundary. The reflection coefficient in- 
creases significantly when fa deviates from zero. Several 
improvements on the characteristic boundary condition 
have been proposed in the literature [4]. One approach 
is the finite wave model given in [2]. 

In the finite wave model, an assumption is made on 
the out-going wave angle and, accordingly, the solution 
at the boundary is constructed from the characteristic 
decomposition associated with that particular propaga- 
tion direction. Specifically, the edge integral along the 
non-reflecting boundary in the semi-discrete equation 
(18) (with a x = (j y = 0) for elements with index n = 0, 

<M x )[A^ 2) u£’ m + AL 2) uJ l ’ m ]ds, (41) 

J 1 0 , m 

is replaced by 

<M x )[A^ 2) u^ m ]ds (42) 

( 2 ) 

where, again, A„ is the boundary normal flux Jacobian 
matrix given in (12) while u^ m is the out-going part of 
the solution along the boundary which is constructed as 
follows. Suppose the presumed out-going wave angle is 
9 , defined in the same fashion as <pi in Figure 7. Then 
the flux in the direction of 9 is (cos 0, sin 0) - F with its 
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Jacobian matrix being A g = cos^A* + sin#A 2 - Let Eg 
be a matrix whose columns are eigenvectors of A$, i-e., 

A, =E,AE- 1 , (43) 

where A is the diagonal eigenvalue matrix. Then 

u% ra - E,I+E„- l u“’ m (44) 

is the solution to be used in the edge integral. Here 
I + is a diagonal matrix whose diagonal values are 1 
and 0 according to the corresponding eigenvalues in A 
being positive and negative respectively. This is equiv- 
alent to decomposing the numerical solution u° h ' m into 
eigenvector components of A# and keeping only the out- 
going components in u^j m . In terms of the expansion 
coefficient vector C" m , the finite wave model leads to 
the following equation to determine the reflection coef- 
ficients: 

(B ( 2 ) ®A( | 2 ) )C 0 ’ m + (B' ( 2 ) ®AL 2 ) )C 1 ' m = (B ( 2 ) ®Ai. 2 ) )C% m 

(45) 

where 


C°ji m = (I <g> Efll+E-^C 0 ’" 1 . (46) 



Figure 10. Reflection coefficient of the finite wave model 
boundary condition. Circle: M = 0; square: M = 0.5. 
uih = 0.5. Open symbols are the results of discrete anal- 
ysis and closed symbols are results from direct numerical 
simulation. Dashed line is the prediction according to 
the solution of partial differential equations. 


In figure 10, we show the reflection coefficient of the 
finite wave model boundary condition with respect to 


an out-going acoustic wave. The non-dimensional fre- 
quency is uh = 0.5. The presumed out-going wave angle 
was taken to be 6 = 7r/4 while the actual out-going wave 
angle fa was varied. As we can see, the reflections are 
effectively annihilated for actual out-going waves with 
an angle fa close to 6. 

Conclusions 

We have presented a discrete analysis of three non- 
reflecting boundary conditions for the discontinuous 
Galerkin method. For the Perfectly Matched Layer ab- 
sorbing boundary condition, we computed the discrete 
wavenumbers and eigenfunctions of the semi-discrete 
equation. Our analysis showed that the reflections at 
the interface of Euler and PML domains caused by the 
discretization process were reduced and approached to 
zero as the resolution of the scheme was increased. For 
the characteristic and finite wave model non-reflecting 
boundary conditions, the discrete reflection errors con- 
verged to those predicted by the non-discrete partial 
differential equations and agreed well with results of 
discrete numerical simulations. 
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Al. Definition of matrices in (24) 

The mass and “stiffness” matrices for equation (24) are 
defined as follows: 


Q = { Q( 1 '}lxl , qt't = / pept'dx 
Jq 


Qx = 

= {<tti'}LxL, 

( Wi = 

Jq dx 

Qy- 

= {QU'}LxL, 

qt't = 

[ 

Jq dy 

B (i) 

II 

o- 

Os 

'-v-' 

X 

bft = 

- / (pt.pt' ds 



J r<>) 

B '(i) 

- {9u'}lxl 

< b'w - 

= / p e pj,ds 
J r<o 


where pj, in matrix B' denotes the basis function of the 
exterior element along edge F^ . 

A2. Non-discrete analysis of the characteristic 
and the finite wave model non-reflecting bound- 
ary conditions 

(1) Characteristic non-reflecting boundary boundary 
condition 

In the non-discrete limit, the characteristic boundary 
condition can be shown to be equivalent to applying 
A~u = 0 at the non-reflecting boundary. Using eigen- 
solutions of the Euler equation (1), the reflection coef- 
ficient of an out-going acoustic wave is found to be 

a 1 - cos 4>i 
Pr — , 

1 — COS <p r 

and that of an out-going vorticity wave is 

a .. sin (pi 
Pr — -i i 

1 — COS (pr 

where pi and p r are defined in Figure 7. It is to be 
noted that the incident and reflected angles are related 
as 


(2) Finite wave model 

Let 9, defined in the same fashion as pi in Figure 7, be 
the presumed angle of the out-going wave used in the 
construction of the solution at the boundary. Then the 
reflection coefficient of an out-going acoustic wave is 

1 - cos (pi - 9) 

Pr ~ 1 - COS (<f>r ~ 9) 

and the reflection coefficient of an out-going vorticity 
wave is 

= sin(</>j - 9) 
r X - COS (p r - 9) ' 


sin p r _ sin pi 
1 + M cos <p r 1 + M COS (pi 


for the out-going acoustic waves and 

sin <p T tan <pi 

1 + M COS (pr M 

for the out-going vorticity waves. 
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